# List of objective functions (logit under various scoring rules)
# Input: None
# Output: List of six objective functions that can be used for model estimation
getof <- function(){
  # Log score
  ls <- function(b, y, x, sel=1){
    s <- x %*% b
    l <- y*log(logit(s))+((1-y)*log(1-logit(s)))
    if (sel==1) -sum(l) else -l
  }
  # Brier score
  br <- function(b,y,x,sel=1) {
    s <- x %*% b
    if (sel==1) 0.5*sum((y-logit(s))^2) else 0.5*(y-logit(s))^2
  }
  # (negative) spherical score
  sp <- function(b,y,x,sel=1){
    p <- logit(x %*% b)
    den <- sqrt(1-2*p+2*p^2)
    num <- y*p+(1-y)*(1-p)
    if (sel==1) -sum(num/den) else -num/den
  }
  # Boosting score 
  bo <- function(b,y,x,sel=1){
    p <- logit(x %*% b)
    if (sel==1) -sum(-y*sqrt((1-p)/p)-(1-y)*sqrt(p/(1-p))) else -(-y*sqrt((1-p)/p)-(1-y)*sqrt(p/(1-p)))
  }
  # As1 score
  rs <- function(b,y,x,sel=1){
    p <- logit(x %*% b)
    if (sel==1) -sum(y*(log(p)-p+1)-(1-y)*p) else -(y*(log(p)-p+1)-(1-y)*p)
  }
  # As2 score
  ds <- function(b,y,x,sel=1){
    p <- logit(x %*% b)
    if (sel==1) -sum(y*(p-1)+(1-y)*(p+log(1-p))) else -(y*(p-1)+(1-y)*(p+log(1-p)))
  }
  return(list(ls, br, sp, bo, rs, ds))
}

# Function for logistic transformation of input "z"
logit <- function(z) 1/(1+exp(-z))

# Simulate data from DGPs used in Section 3 (summarized in Table 2 of the paper) 
# Input: "dgp" -- Index of DGP (1, 2, 3 or 4); "n" -- sample size 
# Output: List with components "x" (draws for regressor), "y" (draws for predictand), 
# "grid" (sensible grid for plotting), "p" (true probability that y = 1, evaluated at grid points), 
# "p.fct" (function for true probability that y = 1), "cdf.x" (cdf of regressor)
sim <- function(dgp = 1, n = 10000){
  # Draw X and p(Y|X)
  if (dgp == 1){
    x <- runif(n, min = -2.5, max = 2.5)
    cdf.x <- function(s) punif(s, min = -2.5, max = 2.5)
    ptrue <- function(s) logit(-0.5*s + 0.2*s^3)
  } else if (dgp == 2){
    x <- runif(n, min = -2.5, max = 2.5)
    cdf.x <- function(s) punif(s, min = -2.5, max = 2.5)
    ptrue <- function(s) logit(-0.5*s + 0.2*s^2)
  } else if (dgp == 3){
    x <- runif(n, min = -1, max = 4)
    cdf.x <- function(s) punif(s, min = -1, max = 4)
    ptrue <- function(s) logit(-0.5*s + 0.2*s^3)
  } else if (dgp == 4){
    x <- runif(n, min = 0, max = 10)
    cdf.x <- function(s) punif(s, min = 0, max = 10)
    ptrue <- function(s) logit(sqrt(s))
  } 
  # Draw Y
  y <- runif(n) < ptrue(x)
  # Make grid for later plotting
  grid <- seq(from = min(x), to = max(x), length.out = 200)
  return(list(x = x, y = y, grid = grid, p = ptrue(grid), p.fct = ptrue, cdf.x = cdf.x))
}

# Function for Diebold-Mariano test, using HAC standard errors
# Input: "ld", sequence of loss differentials. 
# Output: Diebold-Mariano test statistic
dmtest <- function(ld){
  n <- length(ld)
  aux <- data.frame(ld = ld)
  aux <- lm(ld~1,data=aux)
  v <- unname(NeweyWest(aux)*n)
  return(mean(ld)/sqrt(v/n))  
}